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ABSTRACT 


Parabolic equation models solved using the split-step Fourier (SSF) algorithm, such as 
the Monterey Miami Parabolic Equation model, are commonly used to predict 
underwater sound propagation in deep and shallow water environments. Previous studies 
have shown that the SSF algorithm is very accurate in shallow water when there is no 
density discontinuity between the water column and the sediment, but less effective in the 
presence of realistic density discontinuities due to phase errors that accumulate after a 
few kilometers. In this thesis, the standard density-smoothing approach and an alternative 
hybrid split-step/finite difference method are compared. The goal is to decrease the phase 
errors and increase the model’s long-range accuracy. Different depth meshes and range 
step sizes are implemented in the algorithm to find the optimum results for both 
approaches. It is shown that the density-smoothing method provides better results with 
small range step sizes, while the hybrid method is more effective using longer range step 
sizes. However, the smoothing approach provides a more stable convergence of results, 
whereas the hybrid method solution is more sensitive to change in computational grid 
sizes. A more detailed examination of the density smoothing approach suggests good 
accuracy for a few kilometers, while the hybrid method provides improved agreement 
with a benchmark solution at longer ranges. 
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I. INTRODUCTION 


A. BACKGROUND 

The study of the parabolic equation (PE) approximation to the elliptic Helmholtz 
wave equation dates back to mid-1940s, when Leontovich and Fock introduced the PE 
method to the problem of radio-wave propagation in the atmosphere. Since then, 
propagation modeling using the PE has been successfully applied to microwave 
waveguides, optics, laser-beam propagation, plasma physics, seismic propagation and 
underwater acoustics [1]. Tappert was the first to introduce the PE method for underwater 
acoustic propagation in 1974 [2], and the split-step Fourier (SSF) algorithm to solve the 
PE was developed by Hardin and Tappert in 1973 [3], Since then, many articles have 
been published to improve PE propagation modeling in the ocean. According to a survey 
made in 1990, over a period of 15 years, more than 120 articles and technical reports 
related to PE developments in underwater acoustics had been published, as well as many 
more after that time [1]. 

The main advantage of the parabolic equation is modeling underwater sound 
propagation in the ocean at long-range and in range-dependent environments [4]. Several 
methods and algorithms have been introduced to solve the parabolic equation. The most 
common ones are the SSF algorithm, the implicit finite difference (IFD) algorithm, and 
the finite element (FE) method. There are many advantages and disadvantages to these 
techniques. For SSF, the primary advantage is its speed and simplicity [5]. Also, the SSF 
is efficient for long-range, narrow-angle propagation problems in range-dependent media. 
In environments causing strong bottom interactions, however, it loses some of its 
accuracy. IFD and FE methods give more accurate results in wide-angle, bottom¬ 
interacting, and range-independent media. However, they are generally less efficient and 
less stable than the SSF algorithm [1]. 

The standard parabolic approximation originally developed by Tappert was valid 
only for small angles, and its accuracy was degraded by the phase errors. Several 
attempts were made to reduce these errors and support wider angles of propagation. 
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Thompson and Chapman introduced the wide-angle parabolic equation (WAPE) 
approximation, which was based on Feit and Fleck’s operator splitting [6]. This approach 
treats the square-root operator in a way that differs from the standard parabolic equation 
(SPE). It improved the parabolic equation by increasing the propagation angle from 
around 15 degrees to 40 degrees, reducing the phase errors while retaining the usage of 
the efficient SSF algorithm [4], 

As previously stated, the SSF algorithm is known to be less accurate in the 
presence of strong bottom interactions. In order to treat the density contrasts at the 
interfaces between differing media (e.g., water and sediment), a general smoothing 
function is applied to the interface. This smoothing approach is the primary cause of the 
phase errors due to density discontinuity at the interface for long-range propagation. In 
1996, Yevick and Thomson introduced a hybrid split-step/finite difference (SSF/FD) 
technique for modeling acoustic propagation in the presence of non-uniform density 
profdes [7]. This method offers a practical way to adapt existing SSF-based models by 
separating the operator approximations into density-dependent and density-independent 
components rather than relying on smoothing parameters. An overview of this method is 
presented in this thesis with specific numerical examples showing the improvements in 
the solutions. 

B. PROBLEM STATEMENT 

The Monterey-Miami Parabolic Equation (MMPE) model was developed in the 
mid-1990s and since then has been tested for several existing benchmark scenarios. The 
MMPE method computes the underwater acoustic propagation by using the SSF 
algorithm [8], There are various versions of MMPE, including both 2D and 3D versions, 
versions that incorporate rough surfaces, oceanographic internal waves, and a variety of 
other environmental perturbations. In all cases, however, it applies the WAPE 
approximation of Thomson and Chapman [4]. 

Various techniques have been developed to treat the density discontinuity in the 
bottom interaction in order to increase the accuracy of MMPE. In this thesis, our goal is 
to compare the effects of Tappert’s standard density smoothing approach and Yevick and 
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Thomson’s hybrid split-step/finite difference algorithm [7] in the case of the bottom 
interactions for shallow water. We will perform the propagation in a simple Pekeris 
waveguide defined by Paul Hursky [9] and previously investigated by Owens [8] and 
Smith et al. [10]. 

In order to avoid any issues with implementation details in the existing MMPE 
model, a Matlab version of the SSF algorithm has been developed for this thesis based on 
the same operator approximations as the MMPE model. It was also written only for the 
simple, range-independent Pekeris waveguide that motivated this study. 
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II. THEORY 


A. 


PARABOLIC EQUATION DERIVATION 


We start with the 2-D Helmholtz wave equation in cylindrical coordinates, 


]_ 8 _ 

r dr 


dp 

dr 


d 

+ p — 

dz 


1 dp 
p dz 


+ kpi 2 p = 0 , 


(1) 


where k n = — is the reference wave number, n = — 5 — is the acoustic index of 
c 0 c(r,z) 

refraction, c 0 is the reference sound speed, c(r, z) is the spatially varying acoustic sound 

speed, and p(z) is the density that is assumed to change only with depth. All features of 

the environment are presented in c(r, z) and p(z). We define the PE field function V P 

according to 


'P(r, z) = y[k~rp(r, z) e ' v , 
which leads to the fundamental parabolic equation 

PVT/ 

dr 

where the square-root operator Q is defined as 


„ 9 1 d 

Q = , n ~ H— -z- p — 
V kl dz 


f \ d^ 

yP dZ; 


We can also define the field function as 

'f'= 

which satisfies the PE 


( 2 ) 


(3) 


(4) 


(5) 


dr 


= /* b (G-l)'F. 


( 6 ) 
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Here, the square-root operator takes the more compact form of 


2 = 



1 a 2 

kl dz 2 


(7) 


and the effective index of refraction, rt , is defined as 


~2 2 , 

n = n + - 


2 kl 


1 d 2 p 

p dz 2 


1 dp 2 
\P dz , 


( 8 ) 


The original SSF algorithm is based on (8), which is discussed further in the 
section describing the smoothing approach. Either PE allows for a one-way marching 
algorithm solution of the form 

x F(r+ Ar,z) = exp[zk 0 Ar((9-l)] v F(r,z) (9) 


or similar for 'P and Q. We can define the operators as 


2 . 1 d 2 1 dp d 

-, and y = ~— -—— 


s = n 2 -1, p = 


kl dz 2 


k 0 p dz dz 


The original square root-operator then becomes 


2 1 d 


Q = < "-+—P- 


r i d A 


k z 0 dz 


p dz 




^ 1 dp d 1 5’ ^ 

p 1 dz dz p dz 2 


= ] j l + (n 2 -l) + 
= y/l + S + P + 7 


15' 1 dp d 


kl dz~ k~p dz dz 


( 10 ) 


( 11 ) 


which requires further approximation in order to implement the marching algorithm. 

The MMPE model utilizes the Thomson-Chapman WAPE, defined by the 
operator splitting of Feit and Fleck [6], 

g; =v /rT 7 +vr+ 7 /+vr+f-2. 02) 
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Then the propagator in the marching algorithm takes the form 


exp [ik 0 Ar(Q' - 1)] = exp 

id ^-\/l y H" + n Vl £ — 3 j 


« exp 


exp 

iS^l + ju -l) 

exp^/7>(Vl + s -1 j 

= exp 

[is(jT- sy-i) 

exp 

iS^l + ju -l) 

exp[i£(«-l)] 


(13) 


where S=koAr and the non-commuting terms of the operators e, fi, and y are neglected in 
the second line. As strictly noted by Yevick and Thomson [7], the operation order in (13) 
is very important, and the y term should be applied at the end of the range step, 
following application of the fx tenn. 

The operator currently used in MMPE does not contain a y term, since the 
density terms are built into the effective index of refraction. In this case, we simply 
define 

= (14) 


where.? = « 2 -1. The WAPE approximation does not consider the lowest order cross 
tenns between operators, so the lowest order correction becomes 


02 =Q[ 
= Q[ 


\_ 

8 

]_ 

8 


- (f + ju + y) -s 2 - ju 2 
[^(// + y) + (// + y)^ + 


my+ym\- 


(15) 


For completeness, note that the lowest order correction to the smoothing approach would 
be 


Qi - Q\ 

= &- 


]_ 

8 

1_ 

8 


(<? + //) 2 -£ 2 
[£•// + fj,e\. 



(16) 


In the SSF algorithm, e and /i are just scalar operators in the physical (z) and 
wavenumber spaces ( k : ), respectively. 
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If the y terms and cross tenns are added later as corrections, then we can define an 
intennediate step solution as 


V F (r+ A r, z) = FFT 


' 


( 1 - T ^ 



iS 

1 k: A 


exp 

[f k 0 



: FFT [exp [iS(n- 1)) Tfir, z)] 


, (17) 


where we have replaced ^J\ + ju => II—y and Vl Ts => n , consistent with the 

V K 

application in the ^-domain and z-domain. Note that the complete range step in the 
existing MMPE code is obtained by replacing n => rt , where T'in.z) is replaced by 

v F(r,z). 


B. DENSITY-SMOOTHING APPROACH 


The treatment of the density discontinuity in the MMPE model utilizes a density¬ 
smoothing approach at the water-bottom interface. As presented in (8), instead of the 
standard index of refraction, n, we define an effective index of refraction, n , by 


~2 2 , 

n =n + - 


2 kl 


1 d 2 p 
p dz 2 


1 dp' 
\P dz j 


(18) 


The density discontinuity at the water-bottom interface can be defined as 

P(j) = P w + ( P h - P w ) H(z- z h ), 


(19) 


where p w and pt, are the water column and bottom sediment densities, respectively, and 
are assumed to be constants. The function FI serves as an approximation to the Heaviside 
step function and defines a smooth transition according to 


#(z-zj = #«-) = ! 


c 

1 + tanh(-^—) 


( 20 ) 


where £= z-Zb, and L defines the mixing length. 
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The first and second derivatives of the density are defined in the effective index of 
refraction as 


dp_ 

dz 


, ,dH 

\Pb Pw ) ~ 

dz 


(Ph-Pj 
4 L 


sech 




2 L 


( 21 ) 


and 


d 2 p . , d 2 H 

— = (P h ~Pj - 




dz 1 


(Pb-Pj 

41} 


tanh 


^ Isech 2 ^ ^ 


2 L 


2 L 


( 22 ) 


These expressions are then implemented into (18) to define the effective index of 
refraction. 


C. HYBRID SPLIT-STEP / FINITE DIFFERENCE IMPLEMENTATION 

TECHNIQUE 

The hybrid split-step/finite difference PE algorithm for variable density media 
was introduced by Yevick and Thomson in 1996 [7]. Instead of relying on a smoothing 
function to transition the density discontinuity across the interface, they used a hybrid 
technique to specifically solve for the terms containing the density contrast. This 
approach can easily be implemented into existing SSF codes. To accomplish this, they 
separated the effects of the density by splitting the propagator into density-dependent and 
-independent components, as described in the previous section. In this approach, the 
density propagator is solved by using the finite difference technique and a Pade 
expansion [7]. 


The hybrid approach proposes that Q[ requires an additional operation to (17), 
given by 


v F'(r+ Ar,z) = exp id(^j\ + y -1 j 


'P(r+Ar, z). 


(23) 


If we intend to invoke only the Q[ approximation, then this completes the range step 
calculation, and T' T. Here, the additional operation is approximated by a Pade 
expansion defined as 
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1 +—(l + iS)y 

exp| iS(yli+y-l\ *—y-. 

1 + 1 ( 1 -,% 


( 24 ) 


In order to employ a finite difference approach, we then write (23) as 




'P (r+ A r, z) = 


l + i(H-0)y 


'P(r+ Ar,z), 


(25) 


Since direct application of y is problematic, Yevick and Thomson proposed an alternative 
transverse operator: 


p d ( 1 8 ' 
dz\p dz , 


(26) 


Then y = where both //' and // are well defined operators on T. 


1. Finite Difference Technique 

Consider a discretely sampled environment and associated field where index j=nb 
coincides with the water-bottom interface, j<nb is in the water at the shallower depths 
and j>nb is in the bottom at deeper depths. 


Pi j < nb) = p w ; 

P(i = nb) = ^(p w + p b ); 


(27) 


pi i > »6) = p h . 

p d 

The issue here is the evaluation of the term - 

kg dz 


1 fflF 
p dz 


Assuming a smooth, continuous function, a 3-point centered difference scheme to 
evaluate the derivative is defined by 


2 h 


(28) 
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where f ±l is the function evaluated at the index +1, separated at a distance 2 h, to give an 
approximation to the first derivative of the function at the midpoint index 0. Equivalently, 


if we sampled at half index values, then f ( ' - 


1 

Ti 



When the term f ( ' above is discretized and evaluated with this approximation, we 

obtain 

p 8 (\ gy"! p _1 _ \fc^ 

k\ dzyp dz J k\tSz p\dz ) i p x \dz )_ i 

- 2 —— 2 

L 2 2 

—(f'.-'p.)-— (v 

P\ Pj 

2 2 

In a previous analysis by Smith [11], implemented by Owens in 2014 [8], it was assumed 

that since p is discontinuous at the interface, then the factors p± could simply be 

replaced by their appropriate values above or below the index being evaluated. However, 
it has been detennined that for the finite difference approximation to hold, p must be 
treated as a continuous function in the vicinity of index j. In this case, as was defined by 
Thomson [12], 




P\ = \(Px + A>) > P x = -(Po+P-i)- 
2 ^ 2 ^ 


2, Application of Finite Difference Approximations to Field Operators 

Returning to our field operators p and //', we find 


1 8 2 1 d(dV\ 1 

i/T= — r-T= — - «—r - - - 

kl dz 2 kl dz \ dz ) kl^z \ dz )}_ \ dz 

L 2 2 _ 

(^o AzJ (£ 0 AzJ 
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and 


^=4- 

k\ dz 


Po 


( I S V P ' 


Po 


\p dz ) (Az) 2 


Pi P , 



f > 


1 

1 1 

1 

-Vl- 

-+- 

TV — v 

Pl 

Pi P 1 

Pi 

2 

l 3 i) 

2 


(*o Az) 

Combining (31) and (32) withy = //'-//, we obtain 

/'P = (//'-//)¥ 


(32) 


where 


(k 0 Az)~ 



f > 

f 


Po y 

^1- 


Pi 


l 2 ) 

V 


A L + ^ !L _2 
Pi P i 


\ 

f > 


V + 
T 0 T 

Po 1 
p 1 

^ 1 

J 

l 2 ) 



(A: 0 Azj 


(33) 


P + -^ = 


Po 


+ p, 1 
2 2 

Po 


P- = 


(Pl+Po) 

Po 


P -\ ^(Po+P-l) 

p = p + +p- = — + — = 2p 0 
Pi P i 

2 2 


1 1 

- + - 


(Pl+Po) (Po+P-i) 


(34) 


Here, the definition of p matches the definition of the variable p 0 in Yevick and 
Thomson’s paper. However, we reserve p 0 to designate the values of density at the 
central index j=0. 
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In order to find the solution, we can solve the tridiagonal matrix expression 


AY = BY , 


( 35 ) 


where A 




and B 


\ + ^(\ + iS) r 


. A and B are just identity matrices 


except in the vicinity of z ~ z b , where y is not zero. Note that the field vectors T and 'F 

in (25) are now explicitly written as vectors over depth, and T . The solution of (35) 

completes the range step by applying the density operator term in the Q' operator. 


3. Higher Order Correction Implementation 

Yevick and Thomson [7] also proposed a higher order correction term in Q' 2 
given by 


exp|-y[^(// + y) + (// + y)^ + //y + y//]|. 


(36) 


Since // + / = //', this may be approximated by 


r ,y> 1 l-^-id^ju' + ju's + juy + yju) 

ex P \ -y [*/*' + v' e+ w +rfi ]\ =i —^-■ 

^ ' 1 +—iS{s/i' + /i'£ + juY + y/i) 


(37) 


Note that if y = 0, then this higher order correction becomes simply 


exp 




1-- ^—id(^£ju + /us ) 
1 + l^iS(£jU + jus ) 


(38) 


This tenn can be applied to the existing MMPE approach based on the operator 
approximation Q 2 . 


Yevick and Thomson dropped the juy and yju terms since these give third-order 
derivatives and only the lowest order terms are used in the finite difference 
approximation. Instead, they utilized the correction as 
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exp 


% 14 - iS^s/u' + ju's) 


(39) 


This produces the additional correction 
'F(r+ Ar, z) = exp 


16 


—{sju' + ju's) 


v F'(r+ Ar,z), 


(40) 


where we again invoke the matrix approach 

AW = BW'. 


The elements of the matrix A and B are defined by A = 

. Here, s/u'W follows (32) as, 




(41) 

and 


B = 


+A) 


^=7^2 


(K Az ) 


(k 0 Az) 
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III. NUMERICAL RESULTS 


A. MONTEREY-MIAMI PARABOLIC EQUATION GRID SIZES 

The discretization of the ocean environment is necessary and defined by the mesh 
size (A r, A z) [5]. The depth mesh, A z, is used to detennine the desired fast Fourier 
transform (FFT) size, defined by nz. The relation between A z and nz is defined as 
nz=2z max /Az in the SSF algorithm, where z max is the maximum computational depth and 2 
is the factor to account for the image ocean technique common in SSF models. The 
scaling factor, rzfact , is multiplied by the depth mesh to determine the range step size, A r. 
We get a unifonn grid size when rzfact= 1, while the range step size becomes greater than 
the depth grid size when rzfact> 1, or vice versa. 

In this chapter, we first apply the density-smoothing approach and then the hybrid 
SSF/FD method by using different step sizes and depth meshes. In each case, we compare 
our computed transmission loss with a benchmark solution based on the results computed 
from the normal mode model, Couple07, by Richard Evans [13]. The comparison 
between solutions is based on both signal amplitude and the modulation structure in 
range (which depends on the phase of the solution). In many cases, the amplitude is 
found to agree well, but the phase structure shows disagreement with the benchmark 
solutions. 

After the examination of these two methods, we present a comparison between 
the two to show the convergence, sensitivity, and other differences of the density-contrast 
interface. Finally, for completeness, we show the effect of including the higher order 
correction term. 

For all calculations presented in this thesis, the environment was based on a 
scenario originally defined by Hursky [9], This consists of a Pekeris waveguide of depth 
300 m with a water column sound speed of 1500 m/s overlying a semi-infinite sediment 
half-space of sound speed 1700 m/s. Except for the initial test, in which both the water 
and sediment densities equal 1.0 g/cm , the sediment density is defined as 1.5 g/cm . 
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Attenuation in the sediment is defined as 0.5 dB/A or 0.294 dB/m/kHz, while the water 
column is assumed to be lossless. Pressure calculations were made for a point source at 
depth 180 m transmitting a single frequency of 100 Hz. For comparison, a transmission 
loss trace for a receiver depth of 100 m was extracted out to 20 km range. 

B. NO DENSITY DISCONTINUITY CASE 

As has been stated, the SSF algorithm is known to be very accurate in the 
prediction of transmission loss in the absence of the density discontinuities. In this paper, 
we confirmed the consistency of the SSF algorithm with our benchmark Couple07 
solution for such an environment. Figure 1 shows the comparison of these two solutions 
out to 20 km range. The agreement is found to be excellent. 


No Density Discontinuity 


Couple07 



100 


-| -| o _i_i_i_i_i_i_i_i_i_ 

0 2 4 6 8 10 12 14 16 18 20 

Range(km) 

Figure 1. No Density Discontinuity vs. Couple07 
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c. 


DENSITY-SMOOTHING APPROACH 


In the presence of density discontinuities across boundaries, the SSF algorithm 
utilizes the density-smoothing approach, as defined in the previous section. The mixing 
length, L, is a free variable that can be adjusted. In this work, we examined various 
mixing lengths: l/4ko, l/2ko, 1/ko, 3/2ko, 2/ko, where ko is the reference wavenumber. The 
best results were obtained when L was set to 1/ko. We also tried different depth meshes, 
Az, by adjusting nz, and range step sizes, A r, by changing the scaling factor rzfact. We 
acquired an optimum stable, convergent solution with grid size Az=k/15 and rzfact= 5, 
which corresponded to Az=l m with an FFT size of nz=4096, and range step size Ar=5 
m, or A/3. In order to see the effects of the depth mesh and range step size on the results 
using the density-smoothing approach, Figure 2 shows different rzfact values with fixed 
Az=k/15, and Figure 4 shows different Az values with fixed rzfact= 5. Each result is 
compared to our benchmark solution. 
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Figure 2. 


Density-Smoothing Approach, rzfact= 1, 2, 5, 10, 20 


Figure 2 shows the transmission loss comparisons with fixed value of A z =)J 15. 
The scaling factor rzfact is chosen to be 1, 2, 5 in the first panel and 5, 10, 20 in the 
second panel. If we examine the plot, the smoothing approach matches the benchmark 
solution very well for the first several km for all rzfact values considered. Flowever, 
between 4 and 10 km, a phase error appears and becomes significant, especially after 10 
km. It is worth noting that, for Az < /i /10 ■ convergence appears to occur for rzfact= 5, 
with only very minor changes for smaller rzfact values. 


18 

























Figure 3. Convergence for rzfact= 1, 2, 5, 10, 20 


Figure 3 further illustrates the convergence of transmission loss results between 
15 and 20 km. After examining different rzfact values, it is obvious that we cannot see a 
significant change for rzfact< 5. However, when rzfact> 5, the transmission loss amplitude 
changes slowly. For all cases, the phase error remains consistent and is unaffected by 
variations in rzfact. 
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Figure 4. Density-Smoothing Approach, Az=X/2, a/5, a/ 10, a/1 5, a/40 


Figure 4 shows different values of depth mesh, Az, for a fixed value of rzfact= 5. 
The solution converges for Az < /L/10, and Az = 2/15 provides good agreement with the 
benchmark solution for the first several km. The solution begins to degrade at shorter 
ranges for Az > 2/5 , and when Az<a/2, the agreement is poor at all ranges. Even for the 
best solution, however, the phase error begins to appear at ranges beyond about 8 km. 
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Figure 5. Convergence for Az=A/2, A/5, A/10, A/15, AMO 

The expanded view in Figure 5 shows a stable behavior for Az<A/5. For these 
values, the amplitude stays constant, and we observe slight phase changes at ranges 
beyond 15 km. However, smaller values of Az do not improve the phase match with the 
benchmark solution. In fact, for values of Az<A/15, numerical sensitivities begin to 
degrade the solution in the form of small-scale fluctuations in the result. 

D. HYBRID SPLIT-STEP/FINITE DIFFERENCE ALGORITHM 

For the hybrid SSF/FD approach, as with in the density-smoothing approach, we 
applied different range step sizes and depth meshes to find the optimum solution. Using 
the hybrid approach, we found improved agreement with our benchmark solution at 
longer ranges. We acquired the best result by applying a depth mesh Az=A/15=l m, 
corresponding to an FFT size of nz=4096, and rzfact= 10 or Ar= 10 m. 
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Figure 6. Hybrid Split-Step Method, rzfact= 1, 2, 5, 10, 20 


The first panel in Figure 6 shows the results for rzfact= 1, 2, 5, and the second 
panel shows them for rzfact= 5, 10, 20 with a fixed value of Az=A,/15. It can be concluded 
from the plot that when rzfact< 5, then the transmission loss amplitude and phase error 
remain the same. The phase error begins to appear after 8 km but remains small for all 
rzfact values. 
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Figure 7. Convergence for rzfact= 1, 2, 5, 10, 20 


From the expanded view in Figure 7, an interesting observation is noted. First, as 
in the smoothing approach, the solution appears to converge to a stable result for 
rzfact< 5. However, in this case, a better match with the benchmark result is found just 
prior to convergence when rzfact= 10. Larger values of rzfact give much worse results, so 
the improvement for rzfact= 10 is likely just a coincidence. Therefore, the following 
results are computed for rzfact=5. 
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Figure 8. Hybrid Split-Step Method, Az=a/5, a/10, a/1 5, a/40, a/60 


Figure 8 shows the various solutions at a fixed value of rzfact= 5 for A z= a/ 1 5 to 
A z= a/60. It is obvious from the figure that decreasing the Az values does not improve the 
result for rzfact= 5 but introduces significant variability in the transmission loss 
amplitude. The results show the best match with Az=k/15, where agreement with the 
benchmark solution is good beyond 10 km. For Az > A / 5 (results computed but not 
shown), agreement degrades after a few kilometers, with very poor agreement at all 
ranges for Az « A,. 
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Figure 9. Convergence for Az=A/5, a/ 10, a/1 5, A/40, A/60 


Figure 9 displays an expanded view of the last 5 km of the simulation as the Az 
values decreased. The Az values larger than A/15 cause high phase errors and amplitude 
fluctuations, while we see good agreement with Az=A/15. It seems that the hybrid method 
is very sensitive to Az values and show poor convergence in Az. 

E. COMPARISON OF DENSITY-SMOOTHING APPROACH AND HYBRID 
SPLIT STEP ALGORITHM 

The comparison of the density-smoothing and hybrid split approaches gives 
unique results. Figure 10 shows the optimal solution of the two methods compared to 
benchmark solution Couple07. Both methods give the best result when Az=A/15 and 
converge when rzfact= 5. The agreement with the benchmark is consistent out to about 8 
km for the smoothing approach, while this number is almost 18 km for the hybrid 
method. 
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Both the smoothing approach and the hybrid method tend to converge when rzfact 
is around 5. For both methods, the results are stable when rzfact< 5. In either case, it is 
recommended that rzfact should be between 2 and 5 in order to get accurate results. 
Figure 11 illustrates the last 5 km of the solution where hybrid method shows improved 
agreement with the benchmark solution. The transmission loss amplitudes of these two 
methods are almost the same, but the phase error of the smoothing approach is 
significantly higher than the hybrid method at these ranges. 

One of the main differences between these two methods is their sensitivity to Az, 
or nz. It can be observed from Figure 5 and Figure 9 that the hybrid method is more 
sensitive than the smoothing approach since it changes significantly by varying the Az 
values and does not appear to reach a stable, convergent solution. So while the hybrid 
method appears to produce more accurate results, stability of this method appears to 
remain an issue. 
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Figure 10. Smoothing and Hybrid Methods vs. Couple07 with Az=)J 15 and 
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Comparison of Smoothing and Hybrid Methods 



Figure 11. Convergence of Smoothing and Hybrid Methods with Az=)J 15 and 

Ar=5*Az=m 


F. COMPARISON OF HYBRID AND HYBRID HIGHER ORDER 
CORRECTION 

Equations (15) and (36)-(43) define expressions that provide a higher order 
correction to the square root operator. In Figures 12 and 13, we see the comparison of the 
basic hybrid method and the hybrid higher order correction with the benchmark solution. 
For our case, we can conclude that higher order correction does not improve the solution 
any further. 
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Figure 12. Comparison of Hybrid and Hybrid Higher Order Correction with Az=X/15 

and A r= 1 *Az=)J 15 
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Figure 13. Comparison of Hybrid and Hybrid Higher Order Correction with Az=X /15 
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IV. SUMMARY 


The MMPE model, based on the SSF marching algorithm, has been used to 
predict the sound propagation in deep and shallow water environments for many years. 
Many studies have confirmed MMPE’s accuracy in a wide variety of environments. 
However, in the presence of boundary density discontinuities, phase errors have been 
known to appear due to bottom interactions. This degrades the model’s long-range 
accuracy. In this thesis, we compared results of the split-step density-smoothing and 
hybrid split-step/finite difference methods to treat the density discontinuity at the water- 
bottom interface. Of particular interest was the potential reduction in the phase errors 
using the hybrid approach, which could then be easily implemented into existing versions 
for the MMPE with the goal of improving the model’s long range accuracy. 

For the simple environment examined in this thesis, the SSF was found to give 
results consistent with a benchmark solution out to several kilometers, but beyond that, 
the phase errors accumulated and became significant after 10 km. In our examination, we 
obtained the optimum result by defining a depth mesh, Az=A/15. We also found that 
decreasing the depth mesh does not decrease the phase error significantly, and that the 
method tends to converge when Az is between A/10 and A/15. Similarly, the range step 
size plays an important role in determining the best solution. The optimum result for that 
approach is obtained by a range step size Ar=A/3, corresponding to scaling factor, 
rzfact= 5. 

The hybrid split-step/finite difference method resulted in improved solutions 
compared to our benchmark. This approach produced results that closely matched the 
benchmark results computed using Couple07. The best result was obtained with a depth 
mesh Az=A/15 and range step size Ar=2A/3 corresponding to rzfact= 10. However, 
convergence was achieved at rzfact= 5. It is worth noting that both methods produced 
optimal results for the same mesh sizes. 

For both methods, we found that for a fixed depth mesh of Az«A/1 5, the 
solution converged to a reasonably stable state for A r < 5Az = A/3 . However, for a fixed 
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range step size of A r = A/3 , only the split-step smoothing method appeared to converge 
to a stable state for decreasing values of Az. The hybrid method appeared to reach an 
optimal solution at roughly the same depth mesh scale as the split-step smoothing 
method, Az « A/15 , but appeared to degrade for much smaller values of Az. This suggests 
that the finite difference approach to the density discontinuity is much more sensitive to 
the choice of Az. While the solution is improved, this lack of stability should be treated 
with care. Future work should investigate the differences in sensitivity and convergence 
between these two methods in more detail. 

Finally, we evaluated the influence of the higher order correction in the hybrid 
method. For this simple environment, it did not improve the solution significantly. Due to 
the added complexity of this correction and the lack of significant improvement, its 
inclusion in future updates of the MMPE model is not recommended. 
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APPENDIX. MATLAB CODES FOR HYBRID METHOD 


This Matlab program performs a simple, range-independent calculation of a 
Pekeris waveguide problem as defined by Paul Hursky. The calculation is based on a 
split-step Fourier method which, when using the density-smoothing approach defined by 
Tappert and utilized in MMPE, shows phase errors that accumulate with range. In this 
code, the hybrid method of Yevick and Thomson is implemented in order to test the 
correction of the phase error by eliminating the smoothing approach. 

%% PARAMETERS 
format long 

zbot=300; % bottom interface depth 

D=4*zbot; % max computational depth of real ocean 

zs=180; % source depth 

freq=100; % CW signal 

cw=1500; % sound speed in the ocean 

cb=1700; % sound speed in the bottom 

c0=1500; % reference sound speed 

lambda=cO/freq; % reference wavelength 

dz=lambda/10; % first estimate of dz 

NZ=round(2 * D/dz); % first estimate of NZ 

NZ=nextpow2(NZ); % round up to next power of 2 for FFT size 

NZ=2. A NZ; 

D=(NZ/2)*dz; % max computational depth 
zl=(l/2)*D; % depth to begin spatial domain filtering 
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nb=floor(zbot/dz)+l; % depth index of bottom interface 
rzfact=2; % scaling factor between dr and dz ( 
dr=rzfact*dz; % range step size 

dkz=pi/D; % vertical wavenumber sampling 

kO=(2*pi*freq)/cO; % reference wavenumber 

kl=(3/4)*k0; % wavenumber to begin wavenumber domain filtering 

R0= 1; % reference range 

Rmax=20000; % max computational range 

NR=floor(Rmax/dr); % number of range steps in calculation 

% Density Parameters 

rho_w=1000; % density in water 

rho_b=1500; % density in bottom 

% %Preallocate arrays/matrices 

kz=zeros(l,NZ); 

x=zeros(l,NZ); 

A=zeros(l,NZ); 

B=zeros(l,NZ); 

psi=zeros(l,NZ); 

FK=ones(l,NZ); 

Top=zeros(l,NZ); 

PROPK=zeros( 1 ,NZ); 
z=zeros(l,NZ/2); 

FZ=ones( 1 ,NZ/2); 
cz=zeros(l,NZ/2); 
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rho=zeros( 1 ,NZ/2); 
rho_m=zeros( 1 ,NZ/2); 
rho_t=zeros( 1 ,NZ/2); 
rho_p=zeros( 1 ,NZ/2); 
n_prime_sq=zeros( 1 ,NZ/2); 
atten=zeros( 1 ,NZ/2); 
atten_e=zeros( 1 ,NZ/2); 

Uop=zeros( 1 ,NZ/2); 
PROPZ=zeros(l,NZ/2); 
press=zeros(NZ,NR); 
r=zeros(l,NR); 

TL=zeros(NZ,NR); 

% Hybrid Method Parameters 

se=(kO*dz). A 2; 

delta=kO*dr; 

alpha=( 1/4)*(1-1i* delta); 

beta=( 1/4)*(1+1 i * delta); 

q 1 =(rho_w-rho_b )/(3* rho_w+rho_b); 

q2=(rho_b-rho_w)/(3*rho_b+rho_w); 

q3=(rho_w-rho_b)/(rho_w+3 * rho_b); 

q4=(rho_b-rho_w)/(rho_b+3*rho_w); 

a 1 =( 1 -(alpha/se)* q 1); 

a2=(alpha/se)*ql; 

a3=(alpha/se)*q4; 
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a4=( 1 -(alpha/se) * (q3+q4)); 

a5=(alpha/se)*q3; 

a6=( alpha/se) *q2; 

a7=( 1 -(alpha/se)* q2); 

b 1 =( 1 -(beta/ se) * q 1); 

b2=(beta/se)*ql; 

b3=(beta/se)*q4; 

b4=(l-(beta/se)*(q3+q4)); 

b5=(beta/se)*q3; 

b6=(beta/se)*q2; 

b7=( 1 -(beta/ se) * q2); 

% factors used in algebraic solution of tridiagonal matrix 
factO= 1 ,/(a4-a2 * a3/a 1 -a5 * a6/a7); 
factl=(b3-bl*a3/al)*fact0; 
fact2=(b4-b2*a3/al-b6*a5/a7)*fact0; 
fact3=(b5-b7*a5/a7)*fact0; 

% produce r=0 in kz-domain 

kz=[-NZ/2:(NZ/2-l)]*dkz; 

x=kz/kO; 

A=(l -x. A 2). A (-1/4); 

B=-li*NZ*dkz*sqrt(2/(pi*kO))*sin(kz*zs); 

psi=A.*B; 
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% Wavenumber domain filter 


for ii=l:NZ; 

if -kO<kz(ii) && kz(ii)<-kl ; 

FK(ii)=(cos((2 *pi/kO)* (kz(ii)-k 1))). A 2; 
elseif kl<kz(ii) && kz(ii)<kO; 

FK(ii)=(cos((2 *pi/kO)* (kz(ii)-k 1))). A 2; 
elseif kz(ii)<-kO; 

FK(ii)=0; 
elseif kO<kz(ii); 

FK(ii)=0; 

end 

end 

FK=fftshift(FK); 

psi=fftshift(psi).*FK; % starting field with wavenumber filter applied in FFT 

order 

psi=ifft(psi); % transform to z-domain 
% Diffraction propagator (kz-domain) 

Top= 1 - sqrt( 1 -x. A 2); 

PROPK=exp(-li*kO*dr*Top); 

PROPK=fftshift(PROPK); % put in FFT order for application in SSF algorithm 
% Depth 

z=[0:NZ/2-l]*dz; 
for IZ=l:NZ/2; 
if 0<=z(IZ) && z(IZ)<=zl; 
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FZ(IZ)=1; 

elseif zl<z(IZ) && z(IZ)<=(D-dz); 

FZ(IZ)=(cos((pi/(2 * (D-z 1))) * (z(IZ)-z 1))). A 2; 
end 
end 

% Sound Speed Profile (real ocean part of z-domain) 

cz=size(z); 

cz(l:nb-l)=cw; 

cz(nb)=0.5*(cw+cb); 

cz(nb+l :NZ/2)=cb; 

n_prime_sq=(cO./cz). A 2; % For density smoothing approach use effective index of 
refraction instead 

% Density Profile (real ocean part of z-domain) 

rho=size(z); 

rho( 1 :nb-1 )=rho_w; 

rho(nb)=0.5 *(rho_w+rho_b); 

rho(nb+l :NZ/2)=rho_b; 

% Bottom attenuation 

atten(nb)=0.5*(0.5); % in dB/lambda - applying Heaviside step function 
intermediate value 

atten(nb+l:NZ/2)=0.5; % in dB/lambda 
atte n_c=attcn. * (frcq. /cz)/8.686; % in 1/m 
% Lens propagator (real ocean part of z-domain) 

Uop=( 1 -sqrt(n_prime_sq)); 
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PROPZ=exp(-li*kO*dr*Uop).*exp(-atten_e*dr); 

% Range increament 
for IR=1:NR; 

psi(l:NZ/2)=PROPZ.*psi(l:NZ/2); % lens propagator to real ocean 
psi( 1 :NZ/2)=FZ. *psi( 1 :NZ/2); % z-domain filter to real ocean 
psi(NZ/2+2:NZ)=-psi(NZ/2:-l:2); psi(NZ/2+l)=0; % odd symmetry 
psi=fft(psi); % transfonn to kz-domain - output in FFT order 
psi=psi.*PROPK; % diffraction propagator 
psi=psi.*FK; % kz-domain filter 

psi=ifft(psi); % transform to z-domain - output in FFT order 

% Flybrid method - apply in real ocean only (odd symmetry applied after lens 
propagator) 

if IR==1 

disp('Using hybrid finite difference method for density’) 
end 

tmp=fact 1 *psi(nb-1 )+fact2 *psi(nb)+fact3 *psi(nb+1); 
psi(nb-1 )=b 1/a 1 *psi(nb-1 )+b2/a 1 *psi(nb)-a2/a 1 *tmp; 
psi(nb+1 )=b6/a7 *psi(nb)+b7/a7 *psi(nb+1 )-a6/a7 * tmp; 
psi(nb)=tmp; 

% Compute Pressure Field 
r(lR)=lR*dr; 

press(:, IR )=ps i * sqrtf R0/r( IR)); % store pressure field in FFT order 
end 
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% Transmission Loss 


TL=-20*logl0(abs(press(l:NZ/2,:))); 

zout=100; 

nout=find(z>=zout, 1 ); 
if (zout-z(nout-1 ))<(z(nout)-zout) 
nout=nout-l; 
end 
figure 

plot(r/1000,TL(nout,:),'r');axis ij; 
xlim ([0 20]) 
ylim ([30 110]) 
hold on 

plot(rng,tlz,’k');axis ij; 
xlim ([0 20]) 
ylim ([30 110]) 

legend(’PE-S SF-YT-1 stOrder’,'Couple07') 

title(['dr=’ num2str(dr) ’m, lambda/dz=’ num2str(lambda/dz) c0= num2str(c0) 
Ws, rhob= num2str(rho_b) ’kg/m A 3 ]) 
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